﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "seidel")]
    public static unsafe int seidel(IntPtr a_ptr, IntPtr b_ptr, int n, IntPtr x_ptr, double eps)
    {
        double* a = (double*)a_ptr.ToPointer();
        double* b = (double*)b_ptr.ToPointer();
        double* x = (double*)x_ptr.ToPointer();

        return seidel(a, b, n, x, eps);
    }

    /// <summary>
    /// seidel迭代法
    /// </summary>
    /// <param name="a">a[n][n]系数矩阵</param>
    /// <param name="b">b[n]常数向量</param>
    /// <param name="n"></param>
    /// <param name="x">x[n]返回满足精度要求的解向量。若系数矩阵非对角优势，返回解向量0。</param>
    /// <param name="eps">精度要求</param>
    /// <returns>若系数矩阵非对角优势，则显示错误信息，并返回0标志值。否则返回非0标志值。</returns>
    public static unsafe int seidel(double* a, double* b, int n, double* x, double eps)
    {
        int i, j, u, v;
        double p, t, s, q;

        for (i = 0; i <= n - 1; i++)
        {
            u = i * n + i;
            // 置解向量初值
            p = 0.0;
            x[i] = 0.0;
            for (j = 0; j <= n - 1; j++)
            {
                if (i != j)
                {
                    v = i * n + j;
                    p = p + Math.Abs(a[v]);
                }
            }
            //检查系数矩阵是否对角优势
            if (p >= Math.Abs(a[u]))
            {
                //cout << " 系数矩阵非对角优势！" << endl;
                return 0;
            }
        }
        p = eps + 1.0;
        while (p >= eps)
        {
            p = 0.0;
            for (i = 0; i <= n - 1; i++)
            {
                t = x[i]; s = 0.0;
                for (j = 0; j <= n - 1; j++)
                {
                    if (j != i)
                    {
                        s = s + a[i * n + j] * x[j];
                    }
                }
                x[i] = (b[i] - s) / a[i * n + i];
                q = Math.Abs(x[i] - t) / (1.0 + Math.Abs(x[i]));
                if (q > p) p = q;
            }
        }
        return 1;
    }

    /*
      int main()
      { 
          int i;
          double eps;
          double a[4][4]={ 
              {7.0,2.0,1.0,-2.0},
              {9.0,15.0,3.0,-2.0},
              {-2.0,-2.0,11.0,5.0},
              {1.0,3.0,2.0,13.0}};
          double x[4],b[4]={4.0,7.0,-1.0,0.0};
          eps=0.000001;
          i=seidel(&a[0][0], b, 4, x, eps);
          if (i!=0)
          for (i=0;i<4;i++)
            cout <<"x(" <<i <<")=" <<x[i] <<endl;
          return 0;
      }
    */
}
